﻿/*
    Copyright 2010, 2011 Eric Wong
	contact@optimalcycling.com
	http://www.optimalcycling.com
	Optimal Cycling - Advanced power pacing program for cyclists
	
	This file is part of Optimal Cycling.

    Optimal Cycling is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    Optimal Cycling is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    along with Optimal Cycling.  If not, see <http://www.gnu.org/licenses/>.
*/

using System;

namespace OptimalCycling
{
    #region Alias

    using sizeT = System.UInt32;
    using size_array = System.UInt32;

    // Allow us to easily switch between floats and double values
#if REAL_IS_DOUBLE
    using real = System.Double;
    using real2 = Vector2d;
    using real4 = Vector4d;
    using real8 = Vector8d;
#else
    using real = System.Single;
    using real2 = Vector2;
    using real4 = Vector4;
    using real8 = Vector8;
#endif

    #endregion

    /// <summary>
    /// Calculates the air density values utilizing formulas from 
    /// Shelquist, R. (2010, August 2). An Introduction to Air Density and Density Altitude Calculations. 
    /// Retrieved December 30, 2009, from Shelquist Engineering: http://wahiduddin.net/calc/density_altitude.htm
    /// </summary>    
    public static class AirDensity
    {
        // Constants used to calculate air density
        const double P_o = 101325;           // sea level standard pressure, Pa
        const double T_o = 288.15;           // sea level standard temperature, deg K
        const double g = 9.80665;            // gravitational constant, m/s^2
        const double L = 6.5;               // temperature lapse rate, deg K/km
        const double R = 8.31432;            // gas constant, J/mol*deg K
        const double Rd = 287.05;            // for dry air, J/(kg*degK) = 287.05
        const double Rv = 461.495;           // for water vapor, J/(kg*degK) = 461.495
        const double M = 28.9644;            // molecular weight of dry air, gm/mol or kg/kmol
        const double E = 6356000;            // the radius of the earth (for 1976 ISA)

        /// <summary>
        /// Calculates the air density values in kg/m^3 from temperature, pressure, altitude, and dew point temperature.
        /// </summary>
        /// <param name="climate">real number array consisting of: 
        ///     Air Density [kg/m^3], 
        ///     Air Temperature [Celsius], 
        ///     Dew Point Temperature [Celsius], 
        ///     Air Pressure [Pascals]
        ///     
        ///     The other arrays given must have the same length as this array.
        /// </param>
        /// <param name="coursePosLinear">An array containing the linear course positions starting with 0 [metres].</param>
        /// <param name="coursePosPoints">A real4 array containing the x, y, z, 0 position points describing the course where y is the vertical [metres].</param>
        /// <returns>The calculated air density values are stored in the X position of "climate".</returns>
        public static void Calculate(real4[] climate, real[] coursePosLinear, real4[] coursePosPoints)
        {
            for (int i = 0; i < climate.Length; i++)
            {
                double linearPosition = coursePosLinear[i];
                double geopotentialElevation = geometricToGeopotentialAltitude(coursePosPoints[i].Y);
                double positionTemperature = climate[i].Y;

                double saturationVaporPressure = waterVaporPressure(positionTemperature);
                double actualVaporPressure = waterVaporPressure(climate[i].Z);
                double absolutePressure = actualPressure(climate[i].W, geopotentialElevation);
                
                climate[i].X = (real)((absolutePressure / (Rd * (positionTemperature + 273.15))) +
                    (actualVaporPressure / (Rv * (positionTemperature + 273.15))));
            }
        }

        /// <summary>
        /// Converts geometric altitude to geopotential altitude
        /// </summary>
        /// <param name="geometricAltitude">Altitude in metres. The one you would find using a measuring tape.</param>
        /// <returns>The geopotential altitude in metres</returns>
        static double geometricToGeopotentialAltitude(double geometricAltitude)
        {
            return (geometricAltitude * E) / (E + geometricAltitude);
        }

        /// <summary>
        /// Calculates the water vapor pressure of the air as a function of temperature
        /// </summary>
        /// <param name="T">The air temperature in degrees Celsius. 
        /// Use ambient temperature to obtain the saturation vapor pressure.
        /// Use the dew point temperature to obtain the actual vapor pressure.</param>
        /// <returns>The water vapor pressure of the air in Pascals</returns>
        static double waterVaporPressure(double T)
        {
            const double eso = 6.1078;
            const double c0 = 0.99999683;
            const double c1 = -0.90826951E-2;
            const double c2 = 0.78736169E-4;
            const double c3 = -0.61117958E-6;
            const double c4 = 0.43884187E-8;
            const double c5 = -0.29883885E-10;
            const double c6 = 0.21874425E-12;
            const double c7 = -0.17892321E-14;
            const double c8 = 0.11112018E-16;
            const double c9 = -0.30994571E-19;

            double p = (c0 + T * (c1 + T * (c2 + T * (c3 + T * (c4 + T * (c5 + T * (c6 + T * (c7 + T * (c8 + T * (c9))))))))));

            return eso / Math.Pow(p, 8);
        }

        /// <summary>
        /// Calculates the actual (station) pressure in Pascals from the altimeterSetting (ie: pressure)
        /// and the geoPotentialElevation
        /// </summary>
        /// <param name="altimeterSetting">Altimeter (pressure) setting in Pascals</param>
        /// <param name="geoPotentialElevation">Geopotential station elevation in metres</param>
        /// <returns>The actual (station) pressure in Pascals</returns>
        static double actualPressure(double altimeterSetting, double geoPotentialElevation)
        {
            const double k1 = 0.190263;
            const double k2 = 8.417286E-5;

            return millibarToPascals(Math.Pow(Math.Pow(pascalsToMillibar(altimeterSetting), k1) - (k2 * geoPotentialElevation), 1 / k1));
        }

        /// <summary>
        /// Converts the pressure in millibars (aka: hPa) to Pascals
        /// </summary>
        /// <param name="millibarPressure">The pressure in millibar</param>
        /// <returns>The pressure in Pascals</returns>
        static double millibarToPascals(double millibarPressure)
        {
            return 100 * millibarPressure;
        }

        /// <summary>
        /// Converts the pressure in Pascals to millibars (aka: hPa)
        /// </summary>
        /// <param name="pascalPressure">The pressure in millibar</param>
        /// <returns>The pressure in millibars</returns>
        static double pascalsToMillibar(double pascalPressure)
        {
            return 0.01 * pascalPressure;
        }

        /// <summary>
        /// Converts the pressure in mm Hg to millibars
        /// </summary>
        /// <param name="mmHgPressure">The pressure in mm Hg</param>
        /// <returns>The pressure in millibars</returns>
        static double mmHgToMillibars(double mmHgPressure)
        {
            return 33.86379 * mmHgPressure;
        }

    }
}